# Alexander F. Gazmararian
# afg2@princeton.edu

# Set seed for replication
set.seed(10)

# Load packages
library(tidyverse)
library(scales)
library(lmtest)
library(sandwich)
library(modelsummary)
library(broom)
library(margins)
library(kableExtra)

# Load functions and themes 
source("code/fun/savefig.r")
source("code/fun/book_theme.r")
source("code/fun/fix_txt.r")

# Load data
dat <- readRDS("data/NatRegQual_Winter21.rds")
dat <- subset(dat, sample == "credibility")
cj <- readRDS("data/NatQual_Conjoint_Winter21.rds")

# Load baseline model specification
source("code/fun/model_spec.r")

# Create Figure 5.2----
level.plot <-
  lm(selected ~ Reversible - 1, cj) %>%
  tidy() %>%
  dplyr::mutate(
    term = dplyr::case_when(
      term == "ReversibleReversible by lawmakers" ~ "Reversible by lawmakers",
      term == "ReversibleNot reversible" ~ "Not reversible",
      term == "ReversibleReversible by leader" ~ "Reversible by leader"
    ),
    term = factor(
      term,
      ordered = TRUE,
      levels = c("Reversible by leader", "Not reversible", "Reversible by lawmakers")
    )
  ) %>%
  ggplot(aes(x=term, y = estimate, label = percent(estimate, accuracy = 1))) +
  geom_col() +
  geom_text(vjust = -.5) +
  book_theme +
  scale_y_continuous(labels = percent, expand = c(0,0), limits = c(0, .6), breaks = seq(0,.6,.2)) +
  labs(x="", y="") +
  theme(plot.margin = margin(t=10))
level.plot
savefig(level.plot, "5.2_figure_conjointlevel", height = 1.75, filepath = "figures/")

# Create Figure 5.3----
# Specify model
f.choice <- selected ~ TaxCreditAmount + TaxCreditDuration + BipartisanSupport + Reversible + AdditionalInvestmentinCountySchools
# Estimate
ols.choice <-
  lm(f.choice, cj) %>%
  coeftest(., vcov. = vcovCL(., ~ ResponseId)) %>%
  tidy(conf.int = TRUE)
ols.choice <- ols.choice %>%
  mutate(
    level = case_when(
      term %in% c("TaxCreditAmount$15 million per year", "TaxCreditAmount$20 million per year") ~ "Credit Amount",
      term %in% c("TaxCreditDuration10 years", "TaxCreditDuration15 years") ~ "Credit Duration",
      term %in% c("BipartisanSupportHigh", "BipartisanSupportLow") ~ "Bipartisan Support",
      term %in% c("ReversibleNot reversible", "ReversibleReversible by leader") ~ "Reversible",
      term %in% c("AdditionalInvestmentinCountySchoolsGuaranteed $15 million") ~ "Additional Investment in County Schools"
    ),
    level = factor(
      level,
      ordered = TRUE,
      levels = c(
        "Reversible",
        "Bipartisan Support",
        "Credit Amount",
        "Credit Duration",
        "Additional Investment in County Schools"
      )
    ),
    term = case_when(
      term == "TaxCreditAmount$15 million per year" ~ "$15 million per year",
      term == "TaxCreditAmount$20 million per year" ~ "$20 million per year",
      term == "BipartisanSupportHigh" ~ "High",
      term == "BipartisanSupportLow" ~ "Low",
      term == "ReversibleNot reversible" ~ "Not reversible",
      term == "ReversibleReversible by leader" ~ "Reversible by leader",
      term == "AdditionalInvestmentinCountySchoolsGuaranteed $15 million" ~ "Guaranteed $15 million",
      term == "TaxCreditDuration10 years" ~ "10 years",
      term == "TaxCreditDuration15 years" ~ "15 years"
    )
  )
ols.choice <- rbind(ols.choice, data.frame(term="Baseline: $10 million per year",estimate=0,std.error=NA,statistic=NA,p.value=NA,conf.low=NA,conf.high=NA,level="Credit Amount"))
ols.choice <- rbind(ols.choice, data.frame(term="Baseline: 5 years",estimate=0,std.error=NA,statistic=NA,p.value=NA,conf.low=NA,conf.high=NA,level="Credit Duration"))
ols.choice <- rbind(ols.choice, data.frame(term="Baseline: Medium",estimate=0,std.error=NA,statistic=NA,p.value=NA,conf.low=NA,conf.high=NA,level="Bipartisan Support"))
ols.choice <- rbind(ols.choice, data.frame(term="Baseline: Reversible by lawmakers",estimate=0,std.error=NA,statistic=NA,p.value=NA,conf.low=NA,conf.high=NA,level="Reversible"))
ols.choice <- rbind(ols.choice, data.frame(term="Baseline: $20 million (75 percent receive funds)",estimate=0,std.error=NA,statistic=NA,p.value=NA,conf.low=NA,conf.high=NA,level="Additional Investment in County Schools"))
ols.choice$lb90 <- with(ols.choice, estimate + qnorm(0.05) * std.error)
ols.choice$ub90 <- with(ols.choice, estimate + qnorm(0.95) * std.error)
ols.choice <- subset(ols.choice, !is.na(term))
ols.choice$term <- factor(ols.choice$term, ordered = TRUE, levels = c(
  "Baseline: $10 million per year",
  "$15 million per year",
  "$20 million per year",
  "Baseline: $20 million (75 percent receive funds)",
  "Guaranteed $15 million",
  "Baseline: Medium",
  "Low",
  "High",
  "Baseline: 5 years",
  "10 years",
  "15 years",
  "Baseline: Reversible by lawmakers",
  "Not reversible",
  "Reversible by leader"
))
choice.plot <-
  ols.choice %>%
  dplyr::filter(term!="(Intercept)") %>%
  ggplot(aes(x=term,y=estimate)) +
  geom_hline(yintercept=0,color="grey",lty="dashed") +
  geom_errorbar(aes(ymin=conf.low,ymax=conf.high),width=0) +
  geom_errorbar(aes(ymin=lb90,ymax=ub90),width=0,size=1.75) +
  geom_point(size=4) +
  facet_wrap(~level,nrow=5,scale="free_y") +
  coord_flip() +
  labs(x="",y="Change in probability of support",
       caption = "Notes: Thick and thin bars denote 90 and 95 percent confidence intervals. Standard errors clustered by respondent.") +
  scale_x_discrete(label=scales::label_wrap(30)) +
  scale_y_continuous(limits=c(-.175,.175), expand = c(0,0), breaks = seq(-.15,.15,.05),
                     label=scales::percent_format(accuracy=1)) +
  book_theme
choice.plot
savefig(choice.plot, "5.3_figure_conjointchoice", height = 4.8, filepath = "figures/")

# Compare High and Low Credibility Policies----
# Create function for bootstrapped standard errors
bootstrap_high_vs_low <- function(i) {
  # Sample n-1 respondents
  df <- cj %>% dplyr::sample_n(11929, replace = TRUE)
  # Estimate model
  choice.ols <- lm(f.choice, df)
  # Create prediction data frame for the high and low credibility choices
  predict.data <-
    data.frame(
      TaxCreditAmount = c("$10 million per year", "$20 million per year"),
      TaxCreditDuration = c("5 years", "15 years"),
      BipartisanSupport = c("Low", "High"),
      Reversible = c("Reversible by lawmakers", "Not reversible"),
      AdditionalInvestmentinCountySchools = c("$20 million (75% receive funds)", "Guaranteed $15 million")
    )
  # Estimate predicted values
  choice.pred <- predict(choice.ols, predict.data, type="response")
  # Extract predicted probability of selection for the high and low credibility choices
  high.cred <- choice.pred[2]
  low.cred <- choice.pred[1]
  return(data.frame(high.cred, low.cred))
}
# Loop over the function for 1000 bootstrapped samples with replacement
cred_delta <- lapply(1:1000, bootstrap_high_vs_low)
cred_delta_df <- do.call("rbind", cred_delta)
cred_delta_df$id <- 1:nrow(cred_delta_df)
cred_delta_df <-
  cred_delta_df %>%
  tidyr::pivot_longer(cols = c(high.cred, low.cred)) %>%
  dplyr::group_by(name) %>%
  dplyr::summarise(
    lb90 = quantile(value, 0.05),
    ub90 = quantile(value, 0.95),
    lb95 = quantile(value, .025),
    ub95 = quantile(value, .975),
    mean = mean(value)
  )
cred_delta_df

# Heterogeneous effects by partisanship for online appendix----
# Specify Variables
varlist <- c("TaxCreditAmount", "TaxCreditDuration", "BipartisanSupport", "Reversible", "AdditionalInvestmentinCountySchools")
# Estimate model
acie.pid <-
  lm(update(f.choice, ~ . * PartySummary), cj) %>%
  margins_summary(., variables = varlist, at=list(PartySummary=c("Democrat","Neither","Republican")), vcov = vcovCL(., ~ ResponseId)) %>%
  data.frame()
acie.pid
acie.pid.df <-
  acie.pid %>%
  dplyr::mutate(
    attribute = dplyr::case_when(
      factor %in% c("TaxCreditAmount$15 million per year", "TaxCreditAmount$20 million per year") ~ "Credit Amount",
      factor %in% c("TaxCreditDuration10 years", "TaxCreditDuration15 years") ~ "Credit Duration",
      factor %in% c("BipartisanSupportHigh", "BipartisanSupportLow") ~ "Bipartisan Support",
      factor %in% c("ReversibleNot reversible", "ReversibleReversible by leader") ~ "Reversible",
      factor %in% c("AdditionalInvestmentinCountySchoolsGuaranteed $15 million") ~ "Additional Investment in County Schools"
    ),
    level = dplyr::case_when(
      factor == "TaxCreditAmount$15 million per year" ~ "$15 million per year",
      factor == "TaxCreditAmount$20 million per year" ~ "$20 million per year",
      factor == "BipartisanSupportHigh" ~ "High",
      factor == "BipartisanSupportLow" ~ "Low",
      factor == "ReversibleNot reversible" ~ "Not reversible",
      factor == "ReversibleReversible by leader" ~ "Reversible by leader",
      factor == "AdditionalInvestmentinCountySchoolsGuaranteed $15 million" ~ "Guaranteed $15 million",
      factor == "TaxCreditDuration10 years" ~ "10 years",
      factor == "TaxCreditDuration15 years" ~ "15 years"
    )
  )
acie.pid.df$factor <- NULL
acie.pid.df <- rbind(acie.pid.df, data.frame(PartySummary=NA,AME=0,SE=NA,z=NA,p=NA,lower=NA,upper=NA,attribute="Credit Amount",level="Baseline: $10 million per year"))
acie.pid.df <- rbind(acie.pid.df, data.frame(PartySummary=NA,AME=0,SE=NA,z=NA,p=NA,lower=NA,upper=NA,attribute="Credit Duration",level="Baseline: 5 years"))
acie.pid.df <- rbind(acie.pid.df, data.frame(PartySummary=NA,AME=0,SE=NA,z=NA,p=NA,lower=NA,upper=NA,attribute="Bipartisan Support",level="Baseline: Medium"))
acie.pid.df <- rbind(acie.pid.df, data.frame(PartySummary=NA,AME=0,SE=NA,z=NA,p=NA,lower=NA,upper=NA,attribute="Reversible",level="Baseline: Reversible by lawmakers"))
acie.pid.df <- rbind(acie.pid.df, data.frame(PartySummary=NA,AME=0,SE=NA,z=NA,p=NA,lower=NA,upper=NA,attribute="Additional Investment in County Schools",level="Baseline: $20 million (75 percent receive funds)"))
acie.pid.df$lb90 <- with(acie.pid.df, AME + qnorm(0.05) * SE)
acie.pid.df$ub90 <- with(acie.pid.df, AME + qnorm(0.95) * SE)
acie.pid.df$level <- factor(acie.pid.df$level, ordered = TRUE, levels = c(
  "Baseline: $10 million per year",
  "$15 million per year",
  "$20 million per year",
  "Baseline: $20 million (75 percent receive funds)",
  "Guaranteed $15 million",
  "Baseline: Medium",
  "Low",
  "High",
  "Baseline: 5 years",
  "10 years",
  "15 years",
  "Baseline: Reversible by lawmakers",
  "Not reversible",
  "Reversible by leader"
))
choice.pid.plot <-
  acie.pid.df %>%
  dplyr::filter(level!="(Intercept)" & !is.na(PartySummary)) %>%
  ggplot(aes(x=level,y=AME,color=PartySummary,shape=PartySummary)) +
  geom_hline(yintercept=0,color="grey",lty="dashed") +
  geom_errorbar(aes(ymin=lower,ymax=upper),width=0,position=position_dodge(.6)) +
  geom_errorbar(aes(ymin=lb90,ymax=ub90),width=0,size=1.75,position=position_dodge(.6)) +
  geom_point(size=4,position=position_dodge(.6)) +
  facet_wrap(~attribute,nrow=5,scale="free_y") +
  coord_flip() +
  labs(x="",y="Change in probability of support", shape="") +
  scale_x_discrete(label=scales::label_wrap(30)) +
  scale_y_continuous(limits=c(-.175,.175), expand = c(0,0), breaks = seq(-.15,.15,.05),
                     label=scales::percent_format(accuracy=1)) +
  scale_color_manual(values=c("Blue","Black","Red"),name="") +
  book_theme +
  theme(legend.position="bottom")
choice.pid.plot
savefig(choice.pid.plot, filename = "si_conjoint_pid", filepath = "figures/ch5/", height = 6)

# Interpretation checks----
# Create pretty coefficient names
coefnames <- c(
  "(Intercept)" = "Intercept",
  "reversereversible by lawmakers"="Treatment: Reversible by Lawmakers",
  "reversereversible by the Governor"="Treatment: Reversible by the Governor",
  "reversereversible by the President"="Treatment: Reversible by the President",
  "bipart"="Treatment: Low Bipartisanship",
  "schools"="Treatment: Costly Signal",
  "credit10 years"="Treatment: 10 year credit",
  "credit15 years"="Treatment: 15 year credit",
  "age"="Age",
  " I(age^2)"="Age$^2$",
  "Female"="Woman",
  "CollegeDegree"="College (=1)",
  "Black"="Black",
  "Hispanic"="Hispanic",
  "PartySummaryNeither"="Neither party",
  "PartySummaryRepublican"="Republican",
  "pid7_z_cred" = "Partisan strength",
  "ruralhrsa"="Rural",
  "employfull"="Employed (=1)",
  "schools:pid7_z_cred" = "Costly signal x Partisan strength",
  "schools:Dem"="Costly signal x Democrat",
  "bipart:PartySummaryNeither"="Low bipartisanship x Neither party",
  "bipart:PartySummaryRepublican"="Low bipartisanship x Republican",
  "credit15 years:PartySummaryNeither"="15 years x Neither party",
  "credit15 years:PartySummaryRepublican"="15 years x Republican",
  "credit10 years:PartySummaryNeither"="10 years x Neither party",
  "credit10 years:PartySummaryRepublican"="10 years x Republican",
  "gw_index_cred"="Climate beliefs (index)",
  "green_index_cred"="Green jobs (index)",
  "Income_InclLess than $10,000"="Income: Less than 10,000",
  "Income_Incl$20,000 - $29,999"="Income: 20,000 - 29,999",
  "Income_Incl$30,000 - $39,999"="Income: 30,000 - 39,999",
  "Income_Incl$40,000 - $49,999"="Income: 40,000 - 49,999",
  "Income_Incl$50,000 - $59,999"="Income: 50,000 - 59,999",
  "Income_Incl$60,000 - $69,999"="Income: 60,000 - 69,999",
  "Income_Incl$70,000 - $79,999"="Income: 70,000 - 79,999",
  "Income_Incl$80,000 - $99,999"="Income: 80,000 - 99,999",
  "Income_Incl$100,000 - $119,999"="Income: 100,000 - 119,999",
  "Income_Incl$120,000 - $149,999"="Income: 120,000 - 149,999",
  "Income_Incl$150,000 - $199,999"="Income: 150,000 - 199,999",
  "Income_Incl$200,000 - $249,999"="Income: 200,000 - 249,999",
  "Income_Incl$250,000 - $349,999"="Income: 250,000 - 349,999",
  "Income_Incl$350,000 - $499,999"="Income: 350,000 - 499,999",
  "Income_Incl$500,000 or more"="Income: 500,000 or more",
  "Income_InclPrefer not to say"="Income: Prefer not to say"
)

## Reversibility ---- 

# Estimate models
m <- list()
m[[1]] <- lm(reverse_out ~ reverse, dat)
m[[2]] <- lm(update(f.base, reverse_out ~ reverse + .), dat)
m[[3]] <- lm(update(f.base, reverse_out ~ reverse + gw_index_cred + .), dat)
m[[4]] <- lm(update(f.base, reverse_out ~ reverse + gw_index_cred + green_index_cred + .), dat)
m[[5]] <- lm(update(f.base, reverse_out ~ reverse + gw_index_cred + green_index_cred + ruralhrsa + .), dat)

# Create table
file <- "tables/ch5/ols_conjoint_reversibility.txt"
modelsummary(
  m,
  vcov = "HC2",
  stars = c("*"=.1,"**"=.05,"***"=.01),
  coef_map = coefnames,
  gof_map = c("nobs","adj.r.squared"),
  output = "latex",
  escape = FALSE
) %>%
  cat(., file = file)
fix_txt(file)

## Bipartisanship ----

# Estimate Models
bipart.ols <- list()
bipart.ols[[1]] <- lm(update(f.base, bipart_out_i ~ bipart + . + gw_index_cred), dat)
bipart.ols[[2]] <- lm(update(f.base, bipart_out_i ~ bipart * PartySummary + . + gw_index_cred), dat)
bipart.ols[[3]] <- lm(update(f.base, as.numeric(check_bipart) ~ bipart + . + gw_index_cred), dat)
bipart.ols[[4]] <- lm(update(f.base, as.numeric(check_bipart) ~ bipart * PartySummary + . + gw_index_cred), dat)

# Create table
file <- "tables/ch5/ols_conjoint_bipart.txt"
modelsummary(bipart.ols,
             vcov = "HC2",
             stars = c("*"=.1,"**"=.05,"***"=.01),
             coef_map = coefnames,
             gof_map = c("nobs", "adj.r.squared"),
             escape=FALSE,
             output = "latex"
) %>%
  add_header_above(c(" " = 1, "Binary:" = 2, "Scale:" = 2)) %>%
  cat(., file = file)
fix_txt(file)

## Tax Credit ----

# Set up coefficient name
dat$credit <- paste0(dat$credit, " years")
dat$credit <- relevel(factor(dat$credit), ref = "5 years")
dat$check_credit_scale <- as.numeric(dat$check_credit)

# Estimate linear model
credit.ols <- list()
credit.ols[[1]] <- lm(update(f.base, credit_out_i ~ credit + .), dat)
credit.ols[[2]] <- lm(update(f.base, credit_out_i ~ credit + . + gw_index_cred + green_index_cred), dat)
credit.ols[[3]] <- lm(update(f.base, credit_out_i ~ credit * PartySummary + .), dat)
credit.ols[[4]] <- lm(update(f.base, check_credit_scale ~ credit + .), dat)
credit.ols[[5]] <- lm(update(f.base, check_credit_scale ~ credit + . + gw_index_cred + green_index_cred), dat)
credit.ols[[6]] <- lm(update(f.base, check_credit_scale ~ credit * PartySummary + .), dat)

# Create table
file <- "tables/ch5/ols_conjoint_creditlength.txt"
modelsummary(credit.ols,
             vcov = "HC2",
             stars = c("*"=.1,"**"=.05,"***"=.01),
             coef_map = coefnames,
             gof_map = c("nobs", "adj.r.squared"),
             escape=FALSE,
             output = "latex"
) %>%
  add_header_above(c(" " = 1, "Binary:" = 3, "Scale:" = 3)) %>%
  cat(., file = file)
fix_txt(file)

## Costly Signals ----

# Estimate models
ols.school <- list()
ols.school[[1]] <- lm(check_schools_i ~ age + Female + Black + Hispanic + CollegeDegree + PartySummary + ruralhrsa + schools, dat)
ols.school[[2]] <- lm(check_schools_i ~ age + I(age^2) + Female + Black + Hispanic + CollegeDegree + PartySummary + ruralhrsa + employfull + schools, dat)
ols.school[[3]] <- lm(check_schools_i ~ age + I(age^2) + Female + Black + Hispanic + CollegeDegree + PartySummary + ruralhrsa + employfull + Income_Incl + schools, dat)
ols.school[[4]] <- lm(check_schools_i ~ age + I(age^2) + Female + Black + Hispanic + CollegeDegree + Income_Incl + ruralhrsa + schools, dat)
ols.school[[5]] <- lm(check_schools_i ~ age + I(age^2) + Female + Black + Hispanic + CollegeDegree + ruralhrsa + schools * Dem, dat)

# Create table
file <- "tables/ch5/ols_conjoint_schools.txt"
modelsummary(ols.school,
             vcov = "HC2",
             stars = c("*"=.1,"**"=.05,"***"=.01),
             coef_map = coefnames,
             gof_map = c("nobs", "adj.r.squared"),
             escape=FALSE,
             output = "latex"
) %>%
  cat(., file = file)
fix_txt(file)
